Figure 1 - This figures shows…
library(png)
library(maptools)
Checking rgeos availability: TRUE
library(raster)
library(gam)
Loading required package: splines
Loading required package: foreach
foreach: simple, scalable parallel programming from Revolution Analytics
Use Revolution R for scalability, fault tolerance and more.
http://www.revolutionanalytics.com
Loaded gam 1.14
setwd("~/Desktop/Botero postdoc 2016/Human density and the origins of agriculture/")
par(mar=c(0,0,0,20))
d <- readPNG("Larson_dates.png")
plot(seq(0,18, length.out = 19), seq(0,36, length.out = 19), type="n",ylim=c(0,36),xlim=c(0, 18), xaxt="n")
rasterImage(d, 0,0,18,36, interpolate=TRUE, col=d)
Start_of_early_window <- 16-12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 17-4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 34, 34, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 34, 34, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
These dates are provided in the supplimentary information for the Larson (2014) paper. I’ve copied those values into a .csv table provided here.
domestication_times <- read.csv("Domestication timing larson 2014.csv")
dim(domestication_times)
[1] 77 8
| Region | Species | Start.Exploitation | Finish.Exploitation | Start.predomestication | Finish.predomestication | Start.Domestication | Finish.Domestication |
|---|---|---|---|---|---|---|---|
| Southwest asia | Wheat | 12.00 | 11.25 | 11.25 | 11.00 | 11.00 | 9.00 |
| Southwest asia | Barley | 12.00 | 11.25 | 11.25 | 10.50 | 10.50 | 9.00 |
| Southwest asia | Lentil | 12.00 | 11.00 | 11.00 | 10.50 | 10.50 | 9.00 |
| Southwest asia | Pea | 11.50 | 11.00 | 11.00 | 10.00 | 10.00 | 8.50 |
| Southwest asia | Chickpea | 11.00 | 10.50 | 10.50 | 10.25 | 10.25 | 8.25 |
| Southwest asia | Broadbean | NA | NA | NA | NA | 10.50 | NA |
| Southwest asia | Flax | 12.00 | 9.50 | NA | NA | 9.50 | NA |
| Southwest asia | Olive | 10.00 | 6.00 | NA | NA | 6.00 | NA |
| Southwest asia | Sheep | 12.00 | 10.50 | 10.50 | 9.75 | 9.75 | 8.00 |
| Southwest asia | Goat | 12.00 | 10.50 | 10.50 | 9.75 | 9.75 | 8.00 |
| Southwest asia | Pig | 12.00 | 11.50 | 11.50 | 9.75 | 10.25 | 9.00 |
| Southwest asia | Cattle, taurine | 11.50 | 10.50 | 10.50 | 10.25 | 10.25 | 8.00 |
| Southwest asia | Cat | NA | NA | 10.50 | 4.00 | 4.00 | NA |
| South Asia | Tree cotton | 8.50 | 4.50 | NA | NA | 4.50 | NA |
| South Asia | Rice | 8.00 | 5.00 | 5.00 | 4.00 | 4.00 | 2.50 |
| South Asia | Little millet | NA | NA | NA | NA | 4.50 | NA |
| South Asia | Browntop millet | NA | NA | NA | NA | 4.00 | NA |
| South Asia | Mungbean | NA | NA | 4.50 | 3.50 | 3.50 | 3.00 |
| South Asia | Pigeonpea | NA | NA | NA | NA | 3.50 | NA |
| South Asia | Zebu cattle | 9.00 | 8.00 | NA | NA | 8.00 | 6.50 |
| South Asia | Water buffalo | 6.00 | 4.50 | NA | NA | 4.50 | NA |
| East Asia | Broomcorn Millet | 10.00 | 8.00 | NA | NA | 8.00 | NA |
| East Asia | Foxtail millet | 11.50 | 7.50 | NA | NA | 7.50 | NA |
| East Asia | Rice | 10.00 | 8.00 | 8.00 | 7.50 | 7.50 | 5.00 |
| East Asia | Soybean | 8.50 | 5.50 | NA | NA | 5.50 | 4.00 |
| East Asia | Ramie | NA | NA | NA | NA | 5.25 | NA |
| East Asia | Melon | 7.00 | 4.00 | NA | NA | 4.00 | 3.75 |
| East Asia | Pig | 12.00 | 8.50 | NA | NA | 8.50 | 6.00 |
| East Asia | Silkworm | 7.00 | 5.25 | NA | NA | 5.25 | NA |
| East Asia | Yak | NA | NA | NA | NA | 4.25 | NA |
| East Asia | Horse | 7.50 | 6.75 | 6.75 | 5.50 | 5.50 | 4.00 |
| East Asia | Bactrian Camel | NA | NA | NA | NA | 4.50 | NA |
| East Asia | Duck | 2.50 | 1.00 | NA | NA | 1.00 | NA |
| East Asia | Chicken | 6.00 | 4.00 | NA | NA | 4.00 | NA |
| New Guinea | Banana | 10.00 | 7.00 | 7.00 | 4.00 | 4.00 | NA |
| New Guinea | Taro | 10.00 | 7.00 | 7.00 | 4.00 | NA | NA |
| New Guinea | Yam | 10.00 | 7.00 | 7.00 | 4.00 | NA | NA |
| Africa and Arabia | Date palm | 7.00 | 6.00 | NA | NA | 5.00 | NA |
| Africa and Arabia | Sorghum | 8.00 | 4.00 | NA | NA | 4.00 | NA |
| Africa and Arabia | Pearl millet | NA | NA | NA | NA | 4.50 | 3.50 |
| Africa and Arabia | Fonio | NA | NA | NA | NA | 2.50 | NA |
| Africa and Arabia | Cowpea | NA | NA | NA | NA | 3.75 | NA |
| Africa and Arabia | Hyacinth bean | NA | NA | NA | NA | 3.75 | NA |
| Africa and Arabia | Rice | 3.50 | 2.00 | NA | NA | 2.00 | NA |
| Africa and Arabia | Oil palm | 9.25 | 3.50 | NA | NA | 3.50 | NA |
| Africa and Arabia | Cattle | NA | NA | 9.00 | 7.75 | 7.75 | 6.50 |
| Africa and Arabia | Donkey | 9.00 | 5.50 | NA | NA | 5.50 | 3.50 |
| Africa and Arabia | Dromedary camel | 6.50 | 3.00 | NA | NA | 3.00 | NA |
| Africa and Arabia | Guinea fowl | NA | NA | 2.50 | 1.50 | 1.50 | NA |
| North America | Squash | 6.50 | 5.00 | NA | NA | 5.00 | NA |
| North America | Sunflower | 6.00 | 4.75 | NA | NA | 4.00 | NA |
| North America | Sumpweed | 6.00 | 4.50 | NA | NA | 4.00 | NA |
| North America | Pitseed goosefoot | 4.75 | 3.75 | NA | NA | 3.75 | NA |
| Meso-america | Squash (pepo) | NA | NA | NA | NA | 10.00 | 9.50 |
| Meso-america | Maize | 10.00 | 9.00 | NA | NA | 9.00 | NA |
| Meso-america | Foxtail millet-grass | NA | NA | NA | NA | 6.00 | 4.00 |
| Meso-america | Common bean | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Avocado | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Chile pepper | NA | NA | NA | NA | 3.00 | NA |
| Meso-america | Turkey | NA | NA | NA | NA | 2.00 | NA |
| South America | Chili pepper | NA | NA | NA | NA | 6.00 | NA |
| South America | Peanut | NA | NA | NA | NA | 5.00 | NA |
| South America | Cotton | NA | NA | NA | NA | 6.00 | NA |
| South America | Coca | NA | NA | NA | NA | 8.00 | NA |
| South America | Now-minor root crops (arrowroot, leren) | NA | NA | NA | NA | 9.00 | NA |
| South America | Squash | NA | NA | NA | NA | 10.00 | NA |
| South America | Common bean | NA | NA | NA | NA | 5.00 | NA |
| South America | Lima bean | NA | NA | 8.25 | NA | 6.00 | NA |
| South America | Monioc | NA | NA | NA | NA | 7.00 | NA |
| South America | Sweet potato | NA | NA | NA | NA | 5.00 | NA |
| South America | White potato | 7.00 | 4.50 | NA | NA | 4.00 | NA |
| South America | Quinoa | 5.00 | NA | NA | NA | 3.50 | NA |
| South America | Yam | NA | NA | NA | NA | 5.50 | NA |
| South America | Llama | 10.00 | 6.00 | NA | NA | 6.00 | NA |
| South America | Alpaca | 10.00 | 5.00 | NA | NA | 5.00 | NA |
| South America | Guinea pig | NA | NA | NA | NA | 5.00 | NA |
| South America | Muscovy Duck | NA | NA | NA | NA | 4.00 | NA |
par(mar=c(5,4,6,1))
dates <- unlist(domestication_times[3:8])
hist(dates, breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main="All dates in dataset" )
mtext("This tells us about how evenly our evidence is distributed in time", 3, line=1)
hist(dates, breaks = 22, xlim=c(15,0), xlab="Thousand years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main="All dates in dataset with Larson(2014) date windows")
Start_of_early_window <- 12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
hist(dates, breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.2), border=adjustcolor("cornflowerblue", alpha= 0.9), main="", add=TRUE)
mtext("Early Holocene", 3, line = -1, adj=.3)
mtext("Middle Holocene", 3, line= -1, adj=.6)
par(mfrow=c(2,3), mar=c(4,4,2,0))
dim(domestication_times)
[1] 77 8
specific_dates <- domestication_times[,3:8]
for(i in c(1, 3, 5, 2, 4, 6)){
hist(specific_dates[,i], breaks = 22, xlim=c(15,0), xlab="Thousand years ago", col=adjustcolor("cornflowerblue", alpha= 0.5), border=adjustcolor("cornflowerblue", alpha= 0.9), main= names(specific_dates)[i])
Start_of_early_window <- 12
End_of_early_window_start_of_late_window <- 8.2
End_of_late_window <- 4.2
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(0, 30, 30, 0), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
hist(specific_dates[,i], breaks = 22, xlim=c(15,0), xlab="K years ago", col=adjustcolor("cornflowerblue", alpha= 0.2), border=adjustcolor("cornflowerblue", alpha= 0.9), main="", add=TRUE)
}
I’m creating new rows for this table, combining dates in different ways to make the CDFs below look more authentic. This makes it so that pre-ag always happens before post-ag. What I’ve done is given the later date to the earlier date when those dates are missing.
h <- which(is.na(domestication_times[,3]))
domestication_times <- cbind(domestication_times, rep(NA, length(domestication_times[,1])))
domestication_times[,9] <- domestication_times[,3]
domestication_times[h,9] <- domestication_times[h,7]
colnames(domestication_times)[9] <- "adopt exploitation date"
domestication_times[,10] <- domestication_times[,7]
domestication_times[which(is.na(domestication_times[,10])),10] <- 0
colnames(domestication_times)[10] <- "start of ag"
#save(domestication_times, file="~/Desktop/Human density and the origins of agriculture/Domestication timing larson 2014.Rdata")
I think these are best described by a cummulative distribution, showing how they accumulate over time.
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(levels(domestication_times$Region)[ type_number])
print(match)
print(j)
}
[1] "Africa and Arabia"
[1] 7.00 8.00 4.50 2.50 3.75 3.75 3.50 9.25 7.75 9.00 6.50 1.50
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
[1] "East Asia"
[1] 10.00 11.50 10.00 8.50 5.25 7.00 12.00 7.00 4.25 7.50 4.50 2.50 6.00
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
[1] "Meso-america"
[1] 10 10 6 3 3 3 2
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 4, 7, 8
[1] "New Guinea"
[1] 10 10 10
Empirical CDF
Call: ecdf(maxer - match)
x[1:1] = 0
[1] "North America"
[1] 6.50 6.00 6.00 4.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 0.5, 1.75
[1] "South America"
[1] 6.0 5.0 6.0 8.0 9.0 10.0 5.0 6.0 7.0 5.0 7.0 5.0 5.5 10.0 10.0 5.0 4.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
[1] "South Asia"
[1] 8.5 8.0 4.5 4.0 3.5 3.5 9.0 6.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
[1] "Southwest asia"
[1] 12.0 12.0 12.0 11.5 11.0 10.5 12.0 10.0 12.0 12.0 12.0 11.5 4.0
Empirical CDF
Call: ecdf(maxer - match)
x[1:6] = 0, 0.5, 1, ..., 2, 8
par(mfcol=c(2,5), mar=c(4,0,5,0))
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
#print(j)
plot(0,0, xlim=c(15,0), ylim=c(0,100), ylab="Percent of species that will eventually \n be domesticated in a region", xlab="Thousand years ago", main=levels(domestication_times$Region)[ type_number], type="n", yaxt="n")
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 100 * (c(0, j(seq(0, maxer, length.out=100))))
lines(x_seq, y_seq, ylim=c(-1,1))
polygon(c(0, x_seq), c(0, y_seq), border=adjustcolor("cornflowerblue",alpha=1), col=adjustcolor("cornflowerblue", alpha=0.2))
if(i == 2 | i == 1)axis(2)
if(i == 3)mtext("Cummulative distribution function for the accumulation of domesticates", 3, line=3.8, col="cornflowerblue")
}
par(mfcol=c(2,5), mar=c(4,0,5,0))
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
plot(0,0, type="n", xaxt="n", xlab="", bty="n")
mtext("Percent of species that will eventually \n be domesticated in a region", 2, line=-5, cex=0.5)
for(i in 1:8){
type_number <- i
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
#print(j)
plot(0,0, xlim=c(15,0), ylim=c(0,100), ylab="Percent of species that will eventually \n be domesticated in a region", xlab="Thousand years ago", main=levels(domestication_times$Region)[ type_number], type="n", yaxt="n")
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 100 * (c(0, j(seq(0, maxer, length.out=100))))
lines(x_seq, y_seq, ylim=c(-1,1))
polygon(c(0, x_seq), c(0, y_seq), border=adjustcolor("cornflowerblue",alpha=1), col=adjustcolor("cornflowerblue", alpha=0.2))
abline(v= maxer - quantile(j)[2], col="limegreen", lwd=2)
if(i == 2 | i == 1)axis(2)
if(i == 2)mtext("25%", 3, line=3.5, adj=-1, col="limegreen")
if(i == 3)mtext("Cummulative distribution function for the accumulation of domesticates", 3, line=3.8, col="cornflowerblue")
if(i == 4)mtext("Choose a y to predict an x", 3, line=3.3, col="cornflowerblue")
break_one <- maxer
break_two <- maxer - quantile(j)[2]
polygon(x=c(break_two, break_two, break_one, break_one), y=c(0, 1, 1, 0), col=adjustcolor("cornflowerblue", alpha=0.2), border=adjustcolor("cornflowerblue",alpha=1))
lines(x=c(break_two, break_two), y=c(0,-1), col="cornflowerblue")
abline(h = 25, col="limegreen", lwd=2)
}
Make this a function. There is a choice of two methods here. At the end of this section we need to print the desision we’re passing to the later analyses.
origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2,
2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
as.character(origins$CONTINENT)
[1] "W_African_Sav" "Sudanic_Savan" "West Africa T" "Ethipian plat" NA
[6] "E_North_Ameri" "New_Guinea" "Mesoamerica" "N_Lowland_SA" "NW_Lowland_SA"
[11] "Sava_W_India" "S_India" "Ganges_E_Indi" "Chinese_loess" "Japanese"
[16] "Lower-MiddleY" "South trop ch" NA "Southwes amaz" "C/S_Andes"
#subset_order <- c(1, 2, 3, 5, 6, 8, 9, 10, 11, 12, 17, 18)
subset_order <- c(8, 10, 9, 5, 18, 7, 6, 20, 1, 2, 13, 16)
origins_subset <- origins[subset_order,]
origins_subset$CONTINENT
[1] Mesoamerica NW_Lowland_SA N_Lowland_SA <NA> <NA> New_Guinea
[7] E_North_Ameri C/S_Andes W_African_Sav Sudanic_Savan Ganges_E_Indi Lower-MiddleY
18 Levels: C/S_Andes Chinese_loess E_North_Ameri Ethipian plat Ganges_E_Indi ... West Africa T
origins_subset$name
NULL
library(maps)
map()
map(origins, add=TRUE, fill=TRUE, col=adjustcolor("cornflowerblue", alpha=1))
database does not (uniquely) contain the field 'name'.
map()
d <- readPNG("Larson_origins.png")
rasterImage(d, -180, -90, 180, 110, interpolate=TRUE, col=d)
map(add=TRUE)
map(origins, add=TRUE, fill=TRUE, col=adjustcolor("cornflowerblue", alpha=1))
database does not (uniquely) contain the field 'name'.
# need to reproject
This is obviously a bad projection fit right now.
#subset and reorder origins. This is currently done at the end of the plot but should be moved forward.
# Load data for population density
load("PopD_all_December.rdata")
PopD.ALL
class : RasterStack
dimensions : 288, 720, 207360, 18 (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -60, 84 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
names : fourK, fiveK, sixK, sevenK, eightK, nineK, tenK, elevenK, twelveK, thirteenK, fourteenK, fifteenK, sixteenK, seventeenK, eighteenK, ...
min values : 5.611358e-07, 1.067142e-06, 2.508241e-06, 6.317553e-06, 2.286934e-05, 7.631922e-05, 1.272693e-04, 2.118215e-04, 2.602175e-04, 3.226203e-04, 4.390267e-04, 5.572032e-04, 7.313966e-04, 8.286005e-04, 8.297062e-04, ...
max values : 2.051069, 2.013452, 2.142908, 1.888403, 1.863014, 1.880628, 1.650615, 1.678033, 1.697732, 1.499115, 1.517264, 1.443677, 1.464867, 1.453581, 1.436394, ...
# Extract data to a matrix
Pop <- values(PopD.ALL)
r <- raster(PopD.ALL, 1)
r
class : RasterLayer
dimensions : 288, 720, 207360 (nrow, ncol, ncell)
resolution : 0.5, 0.5 (x, y)
extent : -180, 180, -60, 84 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : fourK
values : 5.611358e-07, 2.051069 (min, max)
We need to justify our decision to use a GAM over other models. This should include citations to back up those arguments.
We should make our decisions very transparent here. We should be able to justify our decision of 3 degrees of freedom over other possible values.
# Read the polygons
origins <- readShapePoly('Origins_updated.shp')
# Extract data
per.origin <- extract(r, origins, cellnumber = TRUE, buffer = 100000)
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
names(per.origin) <- origins@data[, 1]
str(per.origin)
List of 20
$ W_African_Sav: num [1:309, 1:2] 92506 92507 92508 92509 92510 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Sudanic_Savan: num [1:306, 1:2] 99050 99051 99052 99053 99054 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ West Africa T: num [1:427, 1:2] 102609 102610 102611 103326 103327 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Ethipian plat: num [1:275, 1:2] 106281 106282 106283 106998 106999 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ NA : num [1:195, 1:2] 67404 67405 67406 68119 68120 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ E_North_Ameri: num [1:170, 1:2] 64270 64271 64984 64985 64986 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ New_Guinea : num [1:24, 1:2] 127360 127361 127362 128080 128081 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Mesoamerica : num [1:73, 1:2] 93032 93033 93034 93035 93036 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ N_Lowland_SA : num [1:39, 1:2] 110373 110374 110375 111092 111093 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ NW_Lowland_SA: num [1:24, 1:2] 123319 123320 123321 123322 123323 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Sava_W_India : num [1:34, 1:2] 83312 84031 84032 84749 84750 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ S_India : num [1:18, 1:2] 99153 99154 99872 99873 99874 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Ganges_E_Indi: num [1:92, 1:2] 85494 85495 85496 85497 85498 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Chinese_loess: num [1:84, 1:2] 72598 72599 73318 73319 73320 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Japanese : num [1:36, 1:2] 59681 59682 60401 61121 61843 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Lower-MiddleY: num [1:131, 1:2] 72547 72548 72549 73267 73268 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ South trop ch: num [1:178, 1:2] 84818 84819 84820 84821 84822 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ NA : num [1:258, 1:2] 62493 62494 62495 62496 62497 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ Southwes amaz: num [1:194, 1:2] 137758 137759 137760 137761 137762 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
$ C/S_Andes : num [1:165, 1:2] 137736 137737 137738 137739 137740 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "cell" "value"
par(mfrow=c(4,5), mar=c(0,0,0,0))
for(h in 1:20){
#h <- 3
originI <- Pop[per.origin[[h]][, 1], ]
x_values <- matrix(c(4:21), dim(originI)[1], 18, byrow=TRUE)
x_value_vector <- as.vector(x_values)
y_value_vector <- scale(as.vector(originI))
density_trend <- cbind(x_value_vector, y_value_vector)
plot(density_trend, col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(-10,10), xlim=c(21,4))
points(density_trend, cex=0.5, col=adjustcolor("cornflowerblue", 0.5))
}
par(mfrow=c(4,5), mar=c(0,0,0,0))
for(h in 1:20){
#h <- 3
originI <- Pop[per.origin[[h]][, 1], ]
x_values <- matrix(c(4:21), dim(originI)[1], 18, byrow=TRUE)
x_value_vector <- as.vector(x_values)
y_value_vector <- scale(as.vector(originI))
density_trend <- cbind(x_value_vector, y_value_vector)
plot(density_trend, col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(-4,8), xlim=c(21,4))
points(density_trend, cex=0.5, col=adjustcolor("grey", 0.1))
gammer <- loess(density_trend[,2] ~ density_trend[,1] )
summary(gammer)
lines(density_trend[,1] ,predict(gammer), col="cornflowerblue", lwd=2)
}
span
the parameter α which controls the degree of smoothing.
The size of the neighbourhood is controlled by α (set by span or enp.target). For α < 1, the neighbourhood includes proportion α of the points, and these have tricubic weighting (proportional to (1 - (dist/maxdist)3)3). For α > 1, all points are used, with the ‘maximum distance’ assumed to be α^(1/p) times the actual maximum distance for p explanatory variables.
par(mfrow=c(4,5), mar=c(0,0,0,0))
for(h in 1:20){
#h <- 3
originI <- Pop[per.origin[[h]][, 1], ]
x_values <- matrix(c(4:21), dim(originI)[1], 18, byrow=TRUE)
x_value_vector <- as.vector(x_values)
y_value_vector <- scale(as.vector(originI))
density_trend <- cbind(x_value_vector, y_value_vector)
#points(density_trend, cex=0.5)
gammer <- loess(density_trend[,2] ~ density_trend[,1], span = .5)
summary(gammer)
predict_gam <- predict(gammer, se=TRUE)
plot(density_trend, col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(-1.35,1.35), xlim=c(21,4))
polygon(x=c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window), y=c(-2, 2, 2, -2), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x=c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window), y=c(-2, 2, 2, -2), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
abline(h=0, lty=2, col="grey")
lines(density_trend[,1] ,predict_gam$fit, col="cornflowerblue")
lines(density_trend[,1] , predict_gam$fit + predict_gam$se.fit, col="grey", lty=1)
lines(density_trend[,1] , predict_gam$fit - predict_gam$se.fit, col="grey", lty=1)
}
# need to add a global mean, an everything but the origins mean, and a buffer around the origins mean.
# Function standardization
std <- function(x) {
b <- (x - min(x)) / (max(x) - min(x))
return(rev(b))
}
diff_df <- function(h){
# Calculating mean and
global.means <- global.SD <- list()
for (j in 1:length(per.origin)) {
#print(j)
originI <- Pop[per.origin[[j]][, 1], ]
time <- 21:4
originI <- na.exclude(originI)
b <- apply(originI, 1, std)
nJ <- nrow(originI)
predictions <- matrix(nrow = nJ, ncol = length(time))
colnames(predictions) <- as.character(time)
for(i in 1:nJ) {
# Need to show a gradient of these df values.
model <- gam(b[, i] ~ s(time, df = h))
col <- sample(rainbow(100), 1)
predictions[i, ] <- predict(model)
}
global.means[[j]] <- apply(predictions, 2, mean)
global.SD[[j]] <- apply(predictions, 2, sd)
}
names(global.means) <- paste(names(per.origin), "Means")
names(global.SD) <- paste(names(per.origin), "SD")
return(list(global.means, global.SD))
}
# Confirm that the x-axis is oriented correctly
for_3 <- diff_df(1)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
i <- 1
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,20))
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
originI <- Pop[per.origin[[1]][, 1], ]
#colnames(originI)
#plot(21:4, rep(0, length(21:4)), type="n", ylim=c(0,1), xlim=c(21,4), xaxt="n", yaxt="n")
for(j in 1:18){
points(23- c(jitter(rep(4, length(originI[,j])), 4.5) + j), originI[,j], pch=19, cex=.1)
}
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]], type="b", pch=names(global.means[[i]]))
axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[1], 3, line=-2)
means_matrix <- matrix(rep(NA,19*20), 20, 19)
colnames(means_matrix) <- c("origin", rev(seq(4, 21, by=1)))
means_matrix[,1] <- names(global.means)
for(i in 1:20){
means_matrix[i,2:19] <- global.means[[i]]
}
kable(means_matrix, caption= "Mean values")
| origin | 21 | 20 | 19 | 18 | 17 | 16 | 15 | 14 | 13 | 12 | 11 | 10 | 9 | 8 | 7 | 6 | 5 | 4 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| W_African_Sav Means | 0.275750492966357 | 0.237971703906143 | 0.200504339081811 | 0.165047402018049 | 0.135516039740239 | 0.121748343406075 | 0.132429431895417 | 0.17103635569842 | 0.235948736377181 | 0.321617048802923 | 0.418964971526823 | 0.514676226549989 | 0.591121916502731 | 0.628378724452797 | 0.621163452983342 | 0.582047087468811 | 0.532681459452571 | 0.482909221409947 |
| Sudanic_Savan Means | -0.0186911478417769 | -0.00828615401443568 | 0.00413446846041128 | 0.0215161762523414 | 0.0473073567212149 | 0.0847812201109172 | 0.136209309975754 | 0.20228390605606 | 0.282347797217804 | 0.374044632344897 | 0.472296827728429 | 0.568020597195695 | 0.648225465565203 | 0.697539021870325 | 0.713328901689326 | 0.707327373243822 | 0.697887725176853 | 0.69056153763021 |
| West Africa T Means | -0.0145401490853128 | -0.00333603963380315 | 0.00977013345024412 | 0.027586956435481 | 0.0533165943413057 | 0.0895594750479319 | 0.138174778488521 | 0.199879802842878 | 0.274517388507917 | 0.360833593695506 | 0.45568720488448 | 0.552897898880818 | 0.643252356019774 | 0.715741781485123 | 0.768806992435192 | 0.809952043255802 | 0.850519193881187 | 0.893430170662931 |
| Ethipian plat Means | 0.0516408507031156 | 0.0622592531360626 | 0.0750787355293174 | 0.0933555709623621 | 0.12092332890413 | 0.16140265228906 | 0.217329358000213 | 0.28895304336974 | 0.373944304902601 | 0.467450912361356 | 0.562582640478741 | 0.650493007468181 | 0.722074582042064 | 0.768618080463592 | 0.790630347643354 | 0.797506826970228 | 0.803240018411344 | 0.812789660231236 |
| NA Means | 0.72808776413489 | 0.6723269776978 | 0.613165929035697 | 0.547108730705425 | 0.472515896153583 | 0.393590592158964 | 0.319669766181732 | 0.262708374009385 | 0.233548745835161 | 0.239616095881227 | 0.28323109870948 | 0.359870592475047 | 0.454992276384195 | 0.545528134731924 | 0.616545768777106 | 0.664148269449493 | 0.691583655592282 | 0.708950521834131 |
| E_North_Ameri Means | -0.0424144312128647 | -0.0202618739720356 | 0.00475194649595462 | 0.0365915191653316 | 0.079541401163135 | 0.136467523436837 | 0.207873937132646 | 0.291025513809716 | 0.380397583621781 | 0.468564149521066 | 0.547361392244133 | 0.609663875755146 | 0.654835785797328 | 0.689283124023394 | 0.721880145858871 | 0.761679332218606 | 0.814909663949399 | 0.875153810181115 |
| New_Guinea Means | 0.304105470512681 | 0.28132745701113 | 0.25734789630324 | 0.230996250201583 | 0.201793192056862 | 0.171544911503563 | 0.144548456661618 | 0.127184052406642 | 0.126877035836262 | 0.15052865629083 | 0.203062414319868 | 0.285273581529037 | 0.389014353820968 | 0.496874952215665 | 0.597210277530265 | 0.683771044637899 | 0.754313180721177 | 0.816774228406938 |
| Mesoamerica Means | -0.0315889824677606 | -0.00787599693194719 | 0.0184956204628149 | 0.0511692835844435 | 0.0938250236223921 | 0.148344738459587 | 0.214601573425648 | 0.290380258066549 | 0.372334867390379 | 0.456479966860125 | 0.538354702268212 | 0.613001244481024 | 0.676612820283987 | 0.726625668176205 | 0.767020170677251 | 0.80679724213028 | 0.855043623685109 | 0.908351302981932 |
| N_Lowland_SA Means | 0.208792535197897 | 0.197230324856633 | 0.185079121786267 | 0.171696272010441 | 0.156804931838545 | 0.141342946965234 | 0.128176787393512 | 0.122340507893562 | 0.130435008736839 | 0.159548596595071 | 0.215485053059509 | 0.300036154887671 | 0.405861213383308 | 0.51607541749585 | 0.62000422933541 | 0.712583747106496 | 0.792272847070064 | 0.866211062729842 |
| NW_Lowland_SA Means | 0.326576386801442 | 0.312436135756889 | 0.296990439986673 | 0.278723633997093 | 0.256416348833303 | 0.230169366346516 | 0.201945873000531 | 0.175714919736147 | 0.157651360929879 | 0.155352073191265 | 0.176275702224548 | 0.225642677767534 | 0.302153168560349 | 0.396329228756827 | 0.497310353388214 | 0.594432476311419 | 0.67925446713634 | 0.756553171147025 |
| Sava_W_India Means | -0.0237936944828927 | -0.00818464775151594 | 0.00933560457136895 | 0.0313826871867378 | 0.0607897006722834 | 0.0992817835941126 | 0.147568966381951 | 0.205468601119744 | 0.272443196878899 | 0.347413657258897 | 0.428584371379107 | 0.513039004997093 | 0.595722652725526 | 0.669265301648986 | 0.731180489316986 | 0.788223314024364 | 0.850666362502511 | 0.917274849593167 |
| S_India Means | -0.0222028081854339 | -0.0066002592031116 | 0.0109240264309662 | 0.0332350046742644 | 0.0636053304857559 | 0.104248025366543 | 0.156161312176903 | 0.218860716380081 | 0.291071667647164 | 0.371586637639489 | 0.458699358753368 | 0.549398894270441 | 0.638133144838091 | 0.71736832702396 | 0.786036841922908 | 0.849257343160103 | 0.913846434565112 | 0.979674443478567 |
| Ganges_E_Indi Means | -0.0168749067644347 | -0.00665207444708958 | 0.00516238304586039 | 0.0208799296723999 | 0.0431774421109922 | 0.0742138287818073 | 0.115649991879242 | 0.168563999805984 | 0.233555951882637 | 0.310487513335769 | 0.398078302603856 | 0.492768744130255 | 0.588625507090279 | 0.676929905267031 | 0.753772199676962 | 0.823106046514196 | 0.892776064568763 | 0.964026386513837 |
| Chinese_loess Means | -0.00992922030088753 | -0.00685627234272808 | -0.00263644972242058 | 0.00456650135793812 | 0.0171800768542588 | 0.0379487356429578 | 0.0695048935574124 | 0.114680636092733 | 0.176234630227922 | 0.255681643811977 | 0.353044690591995 | 0.465903186889885 | 0.587480671879438 | 0.70545496017147 | 0.811956397518547 | 0.90571144015868 | 0.990020919073498 | 1.06943656618837 |
| Japanese Means | 0.474649163724578 | 0.443432703865345 | 0.410041741204845 | 0.372273786700216 | 0.329121724356353 | 0.283030955690219 | 0.240090867441392 | 0.209371090119366 | 0.201589325186606 | 0.226840139404003 | 0.291350184139259 | 0.39377817553382 | 0.519726499901562 | 0.642890482156585 | 0.747382040461932 | 0.829407324688329 | 0.892195580878682 | 0.945280306169006 |
| Lower-MiddleY Means | -0.00682865849885968 | -0.00060483620796883 | 0.00671174779516591 | 0.0168353119264841 | 0.0318371795131959 | 0.0536523053498361 | 0.0841840911894777 | 0.125268923860985 | 0.179000307645022 | 0.247353016268821 | 0.331229454779227 | 0.42919493710873 | 0.535517759646172 | 0.639543869787761 | 0.736041145928061 | 0.82613958614389 | 0.915048715480776 | 1.00384399368775 |
| South trop ch Means | -0.00259675023606447 | -0.00369446566682659 | -0.00363747351558766 | -0.000390054740314468 | 0.00904019752001561 | 0.0284916546492706 | 0.0620623767539311 | 0.113139640613902 | 0.183967610626003 | 0.275388397726442 | 0.385135385537217 | 0.507968233105058 | 0.634326208132707 | 0.750251125909705 | 0.848539337462573 | 0.930537024584384 | 1.00231564635787 | 1.06963296926003 |
| NA Means | 0.0366402042172722 | 0.0345449445801444 | 0.0333093635431334 | 0.0345629129341452 | 0.0409323268055525 | 0.0561025771005159 | 0.0840836906480487 | 0.12805515000113 | 0.189442050381472 | 0.269535827938253 | 0.368930130834347 | 0.484509427569726 | 0.606456684388843 | 0.720252238756432 | 0.817753922638379 | 0.899780585526795 | 0.971748613294788 | 1.03906371492767 |
| Southwes amaz Means | 0.207835723880993 | 0.197019933600991 | 0.185689737006914 | 0.173364610555289 | 0.159843367441753 | 0.145872982114604 | 0.134154420414965 | 0.129609533908725 | 0.138844389475622 | 0.168981654647338 | 0.225948263430072 | 0.312180448593161 | 0.422378015176177 | 0.542599962121189 | 0.664662002239541 | 0.784132607465062 | 0.897503322879635 | 1.00679099727573 |
| C/S_Andes Means | 0.0291205853249027 | 0.0314816768181023 | 0.0344713731256695 | 0.0390464429788137 | 0.0464560633526933 | 0.0581169360685076 | 0.0760540915093959 | 0.103031413769306 | 0.142906851931112 | 0.200040121010661 | 0.27794292932455 | 0.377058478195417 | 0.491456226460239 | 0.608384004959301 | 0.721451798462705 | 0.828972489197 | 0.930712639008333 | 1.02909296975965 |
#global.SD
SD_matrix <- matrix(rep(NA,19*20), 20, 19)
colnames(SD_matrix) <- c("origin", rev(seq(4, 21, by=1)))
SD_matrix[,1] <- names(global.SD)
for(i in 1:20){
SD_matrix[i,2:19] <- global.SD[[i]]
}
kable(SD_matrix, caption= "SD values")
| origin | 21 | 20 | 19 | 18 | 17 | 16 | 15 | 14 | 13 | 12 | 11 | 10 | 9 | 8 | 7 | 6 | 5 | 4 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| W_African_Sav SD | 0.186556185173684 | 0.159071135711959 | 0.130974576751884 | 0.101946183025733 | 0.0724951714757066 | 0.0464817364491077 | 0.0281882850407625 | 0.0212718909296576 | 0.0230484233864053 | 0.0260000666669227 | 0.0272778672869396 | 0.0272414655989988 | 0.0286145237100166 | 0.0359445310006679 | 0.0509667284519141 | 0.0718957662099432 | 0.0961865477719618 | 0.12156398251189 |
| Sudanic_Savan SD | 0.0182791317887036 | 0.0120186235797875 | 0.00857406807198337 | 0.0119997090375925 | 0.0200097740504662 | 0.0296208384871706 | 0.0390634218166876 | 0.0465027632519908 | 0.0502676895837266 | 0.0490728899403298 | 0.0422563308391379 | 0.0301155921728213 | 0.016251502147543 | 0.0230519738470708 | 0.0490762671854076 | 0.0798060762063573 | 0.111190090221826 | 0.142489654117473 |
| West Africa T SD | 0.0128551006323539 | 0.0100809743903953 | 0.0110855392306364 | 0.0157356177847031 | 0.0225339267643659 | 0.0302798766216314 | 0.0378753478987503 | 0.0441002370122807 | 0.0478360621585905 | 0.0482472418609378 | 0.045098717097735 | 0.0390673319501593 | 0.0318512113726889 | 0.0254461299782974 | 0.0263553160311764 | 0.0393349506349416 | 0.0572833097318681 | 0.0763495674818114 |
| Ethipian plat SD | 0.198659638310769 | 0.170440069868076 | 0.141737372134567 | 0.112259279537512 | 0.0864755397538964 | 0.0791421681652602 | 0.0986717402217165 | 0.129156174657549 | 0.154623472712395 | 0.166247254796878 | 0.159947337518384 | 0.135808529359533 | 0.0994473010536208 | 0.0618393023227622 | 0.0356919120363483 | 0.0377802530877869 | 0.0546872860266305 | 0.0726120006791913 |
| NA SD | 0.352588159054744 | 0.332443095173091 | 0.310979361000686 | 0.28580375866142 | 0.254547237522766 | 0.216346143796596 | 0.171980702025272 | 0.12388521679558 | 0.0763676236214709 | 0.0403854093134876 | 0.0435516738756018 | 0.0673123480707615 | 0.0842474074682211 | 0.089847619869054 | 0.0865158533690929 | 0.0824620916709309 | 0.0858640772961775 | 0.09604764429718 |
| E_North_Ameri SD | 0.00744774460623942 | 0.00688305183133896 | 0.00730077417171542 | 0.00859083907761313 | 0.0105320687647222 | 0.012800810910875 | 0.0151914913730156 | 0.0175920880323135 | 0.0199689872021909 | 0.0222421032252577 | 0.0241465903979098 | 0.0252306084926904 | 0.0249032846153157 | 0.0225138982726739 | 0.0181697023167219 | 0.0136609854949669 | 0.0121135163669742 | 0.014408338225588 |
| New_Guinea SD | 0.0890975299550152 | 0.0842301251797607 | 0.0792178489696074 | 0.0736641677047341 | 0.0671594485979282 | 0.0596008898697251 | 0.0511674521224338 | 0.0422240170397645 | 0.0331812675142621 | 0.0254299321136062 | 0.0225457007202054 | 0.0268237812941903 | 0.034037094432116 | 0.0390786810653157 | 0.0406181672896982 | 0.0401664674846312 | 0.0399927675349423 | 0.0409680658412061 |
| Mesoamerica SD | 0.0257631275045289 | 0.0187178864900728 | 0.0164812107253495 | 0.0223842305919749 | 0.0343001607587918 | 0.0492203013799705 | 0.0648668376274174 | 0.0787357939871162 | 0.0882653199553832 | 0.0913017452354271 | 0.0865857410678889 | 0.0743419528904328 | 0.0572427399139474 | 0.0406792417457969 | 0.030609493113632 | 0.0315030760258884 | 0.0386540722319049 | 0.0479903934364773 |
| N_Lowland_SA SD | 0.111921322285171 | 0.107282514497114 | 0.102185138216149 | 0.0958207981245285 | 0.0873593968043197 | 0.0763936956444485 | 0.0630414035072374 | 0.0480814037523411 | 0.0329504357663719 | 0.0219371014627934 | 0.0252549283205093 | 0.0397679138235974 | 0.0556648469342838 | 0.0682249344844884 | 0.0766047903929501 | 0.0827208213023969 | 0.0895838817981019 | 0.0974704711325054 |
| NW_Lowland_SA SD | 0.150541353866891 | 0.139669442413694 | 0.127930041006379 | 0.114205838438166 | 0.097868424572761 | 0.0797876044778232 | 0.0631128649695468 | 0.0525350681168884 | 0.0499851737421644 | 0.0520614502101073 | 0.0541664339731254 | 0.054016322802996 | 0.0513832224992821 | 0.0467970105936129 | 0.0411427020149635 | 0.0364198792817815 | 0.035220746530383 | 0.0381254553823437 |
| Sava_W_India SD | 0.00620248739705606 | 0.00409602037065882 | 0.00439841375218349 | 0.00750428764031471 | 0.012286775530684 | 0.0183044382262446 | 0.0252698645356455 | 0.0327879406646454 | 0.0403222428411792 | 0.0470625286607487 | 0.0518738127485427 | 0.0537645121764797 | 0.052235535975566 | 0.0476849083553341 | 0.0415125953320544 | 0.0349790611537332 | 0.029039548251698 | 0.0251114044895353 |
| S_India SD | 0.0085115048531172 | 0.00892906518867892 | 0.00967136329810452 | 0.0109786376784981 | 0.0132686066270033 | 0.0165765384328679 | 0.0204467930723038 | 0.0238711849145664 | 0.0256672972233344 | 0.0255092175722659 | 0.0237276529091089 | 0.0208724202116767 | 0.0181457061716282 | 0.0170301004732513 | 0.0168873351736634 | 0.0164690070990293 | 0.0154340912923534 | 0.0142507144708916 |
| Ganges_E_Indi SD | 0.00996215193394363 | 0.00691767159299303 | 0.00791433587724963 | 0.0127874855691734 | 0.0199162980428461 | 0.0287137566212716 | 0.0388527650554929 | 0.0499222314748722 | 0.0614506116687209 | 0.0728233254374634 | 0.0831495621605526 | 0.0908947825786612 | 0.0944845150273855 | 0.0928545508935936 | 0.0864932764809146 | 0.0768625361267711 | 0.065847656098237 | 0.0558662960341165 |
| Chinese_loess SD | 0.00723695180791585 | 0.00653465042035558 | 0.0067229060136461 | 0.00783510017868662 | 0.00968749896229179 | 0.0121031590893252 | 0.015015434915449 | 0.0186404465048965 | 0.0229918308268338 | 0.0276609710250012 | 0.0323906245272664 | 0.0365700364119769 | 0.0388300984431191 | 0.037141136645161 | 0.0311556397194864 | 0.0229315384793928 | 0.0156852378394716 | 0.0147092261255239 |
| Japanese SD | 0.184696194095272 | 0.167622323188328 | 0.150274927059077 | 0.131797262192788 | 0.11137082376895 | 0.0892794740360921 | 0.0676920587877228 | 0.050106916646986 | 0.0386317289017512 | 0.0331345985067326 | 0.0322763726395578 | 0.0340554941652883 | 0.0358963772505133 | 0.0361379405079019 | 0.035343175991101 | 0.0349568839559365 | 0.036562520951131 | 0.0402092257222996 |
| Lower-MiddleY SD | 0.0109221052239635 | 0.00879404210265073 | 0.0082644102269953 | 0.00985094611053262 | 0.0133492834979398 | 0.0182250521629163 | 0.0238393891749069 | 0.0293168212659387 | 0.0336239347179238 | 0.0358993675826794 | 0.0357122875762433 | 0.033336309677034 | 0.0300883638849522 | 0.0274478338873639 | 0.0261092333651014 | 0.0261920192405151 | 0.0270668753710159 | 0.0286050078015116 |
| South trop ch SD | 0.0164074485323459 | 0.013017339693911 | 0.0101269327626496 | 0.00883966436288355 | 0.0113546594405518 | 0.0176935756740366 | 0.0263781784905157 | 0.0362001386983728 | 0.0458282304601001 | 0.0540734242720755 | 0.0596617127431108 | 0.0616233193553812 | 0.0592826758862225 | 0.0525468261222113 | 0.0423612740395727 | 0.0308335216922742 | 0.0220855802940787 | 0.0229498120978104 |
| NA SD | 0.0809467917305336 | 0.0685832123328493 | 0.0557119272621779 | 0.0416794929850676 | 0.0273556018582627 | 0.0213125874557977 | 0.0336477364276112 | 0.0522922853903463 | 0.0688255532152485 | 0.0796532138108896 | 0.0829295203707789 | 0.0782121775868238 | 0.0673634830545416 | 0.0546060927712949 | 0.0437646882775398 | 0.0377464374007277 | 0.0370323136420871 | 0.0390839743471241 |
| Southwes amaz SD | 0.125774207173705 | 0.115361452592516 | 0.1042880086271 | 0.0916253149336918 | 0.0769919993552501 | 0.0614840306437265 | 0.0481765788042826 | 0.0417025596083742 | 0.0432490383179266 | 0.0477994697641187 | 0.0500604213341766 | 0.0477221942386258 | 0.0415134022061204 | 0.0341998323542236 | 0.0287859842127018 | 0.0282643985012614 | 0.0333826815416954 | 0.0417094659365358 |
| C/S_Andes SD | 0.0727773302277296 | 0.06676048268613 | 0.0609678447564098 | 0.055024169996668 | 0.0491091885343267 | 0.0446462966444592 | 0.0438008333720804 | 0.0470055485124615 | 0.0524694588100196 | 0.0582765439649155 | 0.0633094382982524 | 0.0669266582058461 | 0.068773937206034 | 0.0681584149345244 | 0.0650020480377401 | 0.0612611199316851 | 0.059915336729057 | 0.0611532400816508 |
par(mfrow=c(4,5), mar=c(0,0,0,0))
for(j in 1:20){
plot(means_matrix[1,2:19], col=adjustcolor("cornflowerblue", alpha=0.8), pch=names(means_matrix[1,2:19]), type="n", xlab="year", ylab="Density", xaxt="n", ylim=c(0,1))
originI <- Pop[per.origin[[j]][, 1], ]
for(h in 1:18){
points(23- c(jitter(rep(4, length(originI[,h])), 4.5) + h), originI[,h], pch=19, cex=.1)
}
lines(means_matrix[j,2:19], col=adjustcolor("cornflowerblue", alpha=0.8), pch=names(means_matrix[j,2:19]), type="b", xlab="year", ylab="Density", xaxt="n")
}
axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(1)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using one degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(2)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using two degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(3)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using three degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(4)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using four degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(5)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using 5 degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(10)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using 10 degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(12)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using 12 degree of freedom
par(mfrow=c(4,5), mar=c(0,0,0,0))
for_3 <- diff_df(17)
global.means <- for_3[[1]]
global.SD <- for_3[[2]]
for(i in 1:20){
plot(4:21, global.means[[i]], col=adjustcolor("cornflowerblue", alpha=0.8), xlab="year", ylab="Density", xaxt="n", type="n", ylim=c(0,1), xlim=c(0,22), xaxt="n", yaxt="n")
polygon(x=22-c(Start_of_early_window, Start_of_early_window, End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("limegreen", alpha= 0.2), border=adjustcolor("limegreen", alpha= 0.9))
polygon(x= 22-c( End_of_early_window_start_of_late_window, End_of_early_window_start_of_late_window, End_of_late_window, End_of_late_window) , y=c(-1, 2, 2, -1), col=adjustcolor("firebrick", alpha= 0.2), border=adjustcolor("firebrick", alpha= 0.9))
polygon(c(4:21,21:4) -3 , c(global.means[[i]] + abs(global.SD[[i]]), rev(global.means[[i]] - abs(global.SD[[i]]))), col="cornflowerblue")
lines(global.means[[i]])
#axis(1, at=seq(1,18, by=1), label=rev(seq(4, 21, by=1)))
mtext(names(global.means)[i], 3, line=-1, cex=.5, adj=.3)
}
GAM model using 17 degree of freedom
# Load patricks productivity PCA data
load('Productivity_ALL.RDATA')
# Load origin shapefiles
origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2,
2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
# Extract the data
prod.origin <- extract(Productivity.ALL, origins)
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
# Mean and SD per region
means <- lapply(prod.origin, colMeans, na.rm = TRUE)
sds <- lapply(prod.origin, sd, na.rm = TRUE)
names(means) <- origins@data$CONTINENT
ymax <- max(unlist(means))
ymin <- min(unlist(means))
time <- 4:21
# Plot
#pdf("productivity.pdf", 20, 30)
par(mfrow = c(5, 4), mar = c(2, 2, 2, 0))
for (i in 1:length(means)) {
plot(y = means[[i]], x = time, xlim = c(21, 4), ylim = c(ymin, ymax),
main = names(means)[i], cex.main = 1, cex.lab = 1, cex.axis = 1,
ylab = "Productivity (PCA axis)", xlab = "Thousand of years ago (k)",
pch = 20, lwd = 1, type = "l",
col = c("purple", "green")[origin.time.region[i]])
up <- sds[[i]] + means[[i]]
down <- means[[i]] - sds[[i]]
lines(up ~ time, lty = 2)
lines(down ~ time, lty = 2)
}
#dev.off()
a <- layout(matrix(c(
1, 1, 1, 1, 1, 1, 1, 1,
3, 6, 7, 8, 9, 10, 11, 4,
3, 5, 5, 5, 5, 5, 5, 4,
3, 12, 13, 14, 15, 16, 17, 4,
2, 2, 2, 2, 2, 2, 2, 2
), 5, 8, byrow=TRUE), width=c(1, 1, 1, 1, 1, 1, 1, 1), height=c(0.5, 1, 1.5, 1, 0.5))
layout.show(a)
frameplot <- function(){
plot(21:0,rep(0, 22), xlim=c(17,4), ylim=c(0, 2.25), type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
frameplot_bottom <- function(){
plot(21:0,rep(0, 22), xlim=c(17,4), ylim=c(-0.25, 2), type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
frameplot()
frameplot_bottom()
d <- readPNG("earth.png")
png(file=paste("40962.png",sep=""),width=2000,height=1000, bg="transparent")
par(mar=c(0,0,0,0))
plot(seq(-180, 180, length.out = 19), seq(-90, 90, length.out = 19), type="n",xlim=c(-180, 180),ylim=c(-90, 90), xaxt="n")
rasterImage(d, -180, -90, 180, 90, interpolate=TRUE, col=d)
polygon(x=c(-180,-180, 180,180), y=c(-90, 90, 90, -90), col=adjustcolor("white", alpha=0.1))
#rasterImage(d, -13.5, -13.5, 375, 375, interpolate=TRUE, col=d)
plot(origins_subset, add=TRUE, col=adjustcolor("white", alpha=.8), xaxt="n", border="white", lwd=4) #still need to reproject!!!
dev.off()
null device
1
###################
type_number <- 8
complex_figure <- function(type_number, i, means, sds){
if(i < 6) polygon(x=c(12,12,8.2,8.2), y=c(-1,3,3,-1), col=adjustcolor("cornflowerblue", alpha=0.4), border=NA)
if(i > 5) polygon(x=c(8.2,8.2,4.2,4.2), y=c(-1,3,3,-1), col=adjustcolor("limegreen", alpha=0.4), border=NA)
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 9]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(j)
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- -c(0, j(seq(0, maxer, length.out=100)))
#lines(x_seq, y_seq, type="l", ylim=c(-1,1))
#polygon(c(0, x_seq), c(0, y_seq), border="black", col=adjustcolor("cornflowerblue", alpha=0.5))
#abline(v= maxer - quantile(j)[2])
break_one_1 <- maxer
break_two_1 <- maxer - quantile(j)[2]
# polygon(x=c(break_two_1, break_two_1, break_one_1, break_one_1), y=c(0, 1, 1, 0), col=adjustcolor("cornflowerblue", alpha=0.5), border=NA)
match <- domestication_times[ which(domestication_times$Region == levels(domestication_times$Region)[ type_number]), 10]
maxer <- max(match, na.rm=TRUE)
j <- ecdf(maxer-match)
print(j)
x_seq <- rev(c(0,seq(0, maxer, length.out=100)))
y_seq <- 2+c(0, j(seq(0, maxer, length.out=100)))
#lines(x_seq, y_seq)
#polygon(c(0, x_seq), c(2, y_seq), border="black", col=adjustcolor("limegreen", alpha=0.5))
break_one_2 <- maxer
break_two_2 <- maxer - quantile(j)[2]
# polygon(x=c(break_two_2, break_two_2, break_one_2, break_one_2), y=c(1, 2, 2, 1), col=adjustcolor("limegreen", alpha=0.5), border=NA)
#abline(v=11)
type <- 1
if(type == 1){
x <- c(means[[i]] , means[[i]] + abs(sds[[i]]), means[[i]] - abs(sds[[i]]))
scaled <- scale(x , center=FALSE)
meanss <- scaled[1:18]
sdss_plus <- scaled[19:36]
sdss_minus <- scaled[37:54]
#abline(v=10, col="red")
length(scaled)
#lines(4:21, means[[i]] + sds[[i]])
#polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
polygon(x=c(21:4,4:21), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
}
if(type == 2){
x <- c(means[[i]] , means[[i]] + abs(sds[[i]]), means[[i]] - abs(sds[[i]]))
scaled <- x + 1 #scale(x , center=FALSE)
meanss <- scaled[1:18]
sdss_plus <- scaled[19:36]
sdss_minus <- scaled[37:54]
#abline(v=10, col="red")
length(scaled)
#lines(4:21, means[[i]] + sds[[i]])
polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
}
if(type == 3){
x <- c(means[[i]] , means[[i]] + abs(sds[[i]]), means[[i]] - abs(sds[[i]]))
scaled <- x #scale(x , center=FALSE)
meanss <- scaled[1:18]
sdss_plus <- scaled[19:36]
sdss_minus <- scaled[37:54]
#abline(v=10, col="red")
length(scaled)
#lines(4:21, means[[i]] + sds[[i]])
polygon(x=c(4:21, 21:4), y=c(sdss_plus, rev(sdss_minus)), col=adjustcolor("firebrick", alpha=1), border="white")
}
means_long_y <- c(1,1,1,1,1, meanss)
means_long_x <- c(0:4, 4:21)
break_one <- break_one_2
break_two <- break_two_2
# polygon(x=c(break_one, break_one, 22, 22), y=c(1, 2, 2, 1), col=adjustcolor("white", alpha=0.8), border=NA)
# polygon(x=c(break_two, break_two, break_one, break_one), y=c(1, 2, 2, 1), col=adjustcolor("white", alpha=0), border=NA)
# polygon(x=c(-1,-1, break_two , break_two), y=c(1.9, 3.1, 3.1, 1.9), col=adjustcolor("white", alpha=0.8), border=NA)
#abline(v= break_one, col="white")
#abline(v= break_two, col="white")
break_one <- break_one_1
break_two <- break_two_1
# polygon(x=c(break_one, break_one, 22, 22), y=c(0, 1, 1, 0), col=adjustcolor("white", alpha=0.8), border=NA)
# polygon(x=c(break_two, break_two, break_one, break_one), y=c(0, 1, 1, 0), col=adjustcolor("white", alpha=0), border=NA)
# polygon(x=c(-1,-1, break_two , break_two), y=c(-1.1, .1, .1, -1.1), col=adjustcolor("white", alpha=0.8), border=NA)
#abline(v= break_one, col="white")
#abline(v= break_two, col="white")
#lines(x=c(break_one_2, break_one_2), y=c(1,3), col="white")
#lines(x=c(break_one_1, break_one_1), y=c(1,-1), col="white")
#lines(x=c(break_two_2, break_two_2), y=c(1,3), col="white")
#lines(x=c(break_two_1, break_two_1), y=c(1,-1), col="white")
#lines(4:21, meanss)
lines(21:4, meanss)
}
frameplot()
complex_figure(7, 1, global.means, global.SD)
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 3.5, 4, 4.5
axis(1)
axis(2)
quartz(width=8, height=8)
layout(matrix(c(
1, 1, 1, 1, 1, 1, 1, 1,
3, 6, 7, 8, 9, 10, 11, 4,
3, 5, 5, 5, 5, 5, 5, 4,
3, 12, 13, 14, 15, 16, 17, 4,
2, 2, 2, 2, 2, 2, 2, 2
), 5, 8, byrow=TRUE), width=c(1, 1, 1, 1, 1, 1, 1, 1), height=c(0.5, 1, 1.5, 1, 0.5))
par(mar=c(0,0,0,0))
# 1-4 label margins
blankplot <- function(){
plot(0,0, xlim=c(4,21), ylim=c(1, 1.25), bty="n", type="n", xaxt="n", yaxt="n", xlab="", ylab="")
}
blankplot()
blankplot()
blankplot()
blankplot()
origins <- readShapePoly('Origins_updated.shp')
origin.time.region <- c(2, 2, 1, 1, 1, 2, 2, 1, 2, 2,
2, 2, 1, 2, 2, 2, 2, 2, 2, 2) # 1 = early; 2 = middle
as.character(origins$CONTINENT)
[1] "W_African_Sav" "Sudanic_Savan" "West Africa T" "Ethipian plat" NA
[6] "E_North_Ameri" "New_Guinea" "Mesoamerica" "N_Lowland_SA" "NW_Lowland_SA"
[11] "Sava_W_India" "S_India" "Ganges_E_Indi" "Chinese_loess" "Japanese"
[16] "Lower-MiddleY" "South trop ch" NA "Southwes amaz" "C/S_Andes"
#subset_order <- c(1, 2, 3, 5, 6, 8, 9, 10, 11, 12, 17, 18)
subset_order <- c(8, 10, 9, 5, 18, 7, 6, 20, 1, 2, 13, 16)
origins_subset <- origins[subset_order,]
origins_subset$CONTINENT
[1] Mesoamerica NW_Lowland_SA N_Lowland_SA <NA> <NA> New_Guinea
[7] E_North_Ameri C/S_Andes W_African_Sav Sudanic_Savan Ganges_E_Indi Lower-MiddleY
18 Levels: C/S_Andes Chinese_loess E_North_Ameri Ethipian plat Ganges_E_Indi ... West Africa T
d <- readPNG("earth.png")
png(file=paste("40962.png",sep=""),width=2000,height=1000, bg="transparent")
par(mar=c(0,0,0,0))
plot(seq(-180, 180, length.out = 19), seq(-90, 90, length.out = 19), type="n",xlim=c(-180, 180),ylim=c(-90, 90), xaxt="n")
rasterImage(d, -180, -90, 180, 90, interpolate=TRUE, col=d)
polygon(x=c(-180,-180, 180,180), y=c(-90, 90, 90, -90), col=adjustcolor("white", alpha=0.1))
#rasterImage(d, -13.5, -13.5, 375, 375, interpolate=TRUE, col=d)
plot(origins_subset, add=TRUE, col=adjustcolor("white", alpha=.8), xaxt="n", border="white", lwd=4) #still need to reproject!!!
dev.off()
quartz
2
d <- readPNG("40962.png")
dim(d)
[1] 1000 2000 4
par(mar=c(0,0,0,0))
plot(0:360,0:360,type="n",xlim=c(20,360),ylim=c(65,295), yaxt="n", xaxt="n")
rasterImage(d, -28.5, -13.5, 388, 375, interpolate=TRUE, col=d)
axis(2, label=seq(-90, 90, length.out = 19), at=seq(1, 360, length.out = 19), las=1)
mtext("latitude", 2, line=4, at=180)
abline(h=seq(1, 360, length.out = 19), col=adjustcolor("grey10", alpha= 0.4), lwd=1)
abline(h=180, col=adjustcolor("white", alpha= .5), lwd=1)
load('PopD_all_December.rdata')
# Extract the data
prod.origin <- extract(PopD.ALL, origins_subset)
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
library(matrixStats)
matrixStats v0.51.0 (2016-10-08) successfully loaded. See ?matrixStats for help.
# Mean and SD per region
means <- lapply(prod.origin, colMeans, na.rm = TRUE)
sds <- lapply(prod.origin, colSds, na.rm = TRUE)
## new values from Bruno's GAM model (produced in script called Population_Trend_per_y.R)
means <- global.means
sds <- global.gams
names(means) <- origins_subset@data$CONTINENT
ymax <- max(unlist(means))
ymin <- min(unlist(means))
time <- 4:21
#plot(origins)
#means[[1]] +
#sds[[1]]
#scale(as.numeric(means[[1]]), center=FALSE)
name_vector <- as.character(origins_subset@data$CONTINENT)
for(i in 1:12){
if(i > 6){frameplot()}else{frameplot_bottom()}
## customize polygons for each graph
if(i == 1){ #mesoamerica #values from Larson
complex_figure(3, i, means, sds)
}
#########
if(i == 2 ){ #NW lowlands SA #values from Larson
complex_figure(6, i, means, sds)
}
#########
if( i == 3){ #NW lowlands SA #values from Larson
complex_figure(6, i, means, sds)
}
#########
if(i == 4){ #Fertile crescent aka Southwest asia #values from Larson
complex_figure(8, i, means, sds)
}
#########
if(i == 5){ #loess plateau #values from Larson
complex_figure(2, i, means, sds)
}
#########
if(i == 6){ #new guinea #values from Larson
complex_figure(4, i, means, sds)
}
#########
if(i == 7){ #Eastern N.A. #values from Larson
complex_figure(5, i, means, sds)
}
#########
if(i == 8){ #Andes #values from Larson
complex_figure(6, i, means, sds)
}
#########
if(i == 9){ #W. African Sav #values from Larson
complex_figure(1, i, means, sds)
}
#########
if(i == 10){ #Sudanic sav #values from Larson
complex_figure(1, i, means, sds)
}
#########
if(i == 11){ #Ganges #values from Larson
complex_figure(7, i, means, sds)
}
#########
if(i == 12){ #loess #values from Larson
complex_figure(2, i, means, sds)
}
#lines(4:21, means[[i]])
abline(h = 1, col=adjustcolor("forestgreen", alpha=.5), lty=2)
# add axes to some locations
if(i == 1 | i == 7){axis(2, at=seq(0,2, by=0.25), label=seq(0,2, by=0.25), las=1)}
if(i == 6 | i == 12){axis(4, at=seq(0,2, by=0.25), label=seq(0,2, by=0.25), las=1)}
#if(i == 6 | i == 12){axis(4, at=seq(2,3, by=0.25), label=seq(0,1, by=0.25), las=1)
# axis(4, at=seq(-1,0, by=0.25), label=rev(seq(0,1, by=0.25)), las=1)
# }
if(i > 6){axis(1)} else{axis(3)}
# add text
if(i < 7){polygon(x=c(-30, -30, 30, 30), y=c(-0.1, -0.5, -0.5, -0.1), col="black")
mtext(name_vector[i], 1, line=-1.2, col="white", cex=0.5)}
if(i > 6){polygon(x=c(-30, -30, 30, 30), y=c(2.1, 2.5, 2.5, 2.1), col="black")
mtext(name_vector[i], 3, line=-1.2, col="white", cex=0.5)}
# add axis labels
if(i == 1 | i == 7){mtext("scaled density potential", 2, line=4, at=1)}
if(i == 3){mtext("Thousand years before present", 3, line=3.5, at =5)}
if(i == 9){mtext("Thousand years before present", 1, line=3.5, at =5)
}
}
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 4, 7, 8
Empirical CDF
Call: ecdf(maxer - match)
x[1:5] = 0, 1, 4, 7, 8
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 1, 2, ..., 6, 6.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 1, 2, ..., 6, 6.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:6] = 0, 0.5, 1, ..., 2, 8
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 0.5, 0.75, ..., 5, 7
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 0.5, 1, ..., 4.5, 7.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:1] = 0
Empirical CDF
Call: ecdf(maxer - match)
x[1:2] = 0, 4
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 0.5, 1.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:3] = 0, 1, 1.25
Empirical CDF
Call: ecdf(maxer - match)
x[1:8] = 0, 1, 2, ..., 5, 6
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 1, 2, ..., 6, 6.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 2.25, 2.75, ..., 5.75, 6.25
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.25, 1.25, ..., 6.75, 7.75
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 2.25, 2.75, ..., 5.75, 6.25
Empirical CDF
Call: ecdf(maxer - match)
x[1:7] = 0, 0.5, 1, ..., 5, 5.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:4] = 0, 3.5, 4, 4.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:11] = 0, 0.5, 2, ..., 7.75, 9.5
Empirical CDF
Call: ecdf(maxer - match)
x[1:9] = 0, 0.5, 1, ..., 4.5, 7.5
saveToPDF <- function(...) {
d = dev.copy(pdf,...)
dev.off(d)
}
saveToPNG <- function(...) {
d = dev.copy(png,...)
dev.off(d)
}
## Try them out
saveToPDF("my.pdf", height=8,width=8)
quartz
2
saveToPNG("my.png", height=8, width=8, units="in", res=300)
quartz
2
dev.off()
null device
1